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Abstract: A C++ implementation of the D s -dimensional unitarity cut algorithm for the 
numerical evaluation of the virtual contribution to NLO QCD amplitudes is presented. The 
current version includes an arbitrary number of external gluons with gluonic propagators 
in the loop. The building blocks are tree level color-ordered amplitudes with gluons and 
with gluons and two scalars in five dimensions. Numerical stability issues are addressed 
and agreement has been reached with the results published in the literature. 



1. Introduction 



In anticipation of the LHC data the need for accurate calculations of Standard Model 
processes beyond the leading order is unequivocal [1]. An algorithm that is able to perform 
such next to leading order calculations in a relatively generic and thus automatable way has 
been since long a withstanding goal of the community. Such a possibility seems within reach 
now, after the significant progress of the last two years, triggered by the OPP reduction 
method at the integrand level [2], [3], [4]. This approach has been implemented in a 
publically released code, Cut Tools [3] and has been employed to calculate the six-photon 
amplitude [5] and the NLO QCD corrections to trivector boson production [6]. 

Based on the idea of numerical reduction at the integrand level, the authors of [7] 
employed unitarity-cut methods to recover the cut constructible part of the virtual NLO 
amplitude using tree level amplitudes as building blocks. It was later shown [8] that 
the rational part can also be evaluated numerically by applying the same algorithm in 
higher but integer number of dimensions, D s . The generalization to fermions has been 
demonstrated in [9] and the first phenomenological applications using unitarity cuts have 
been achieved in [10]. 

In parallel, a semi-anaytic approach is engineered by the Black Hat collaboration [11], 
[12], [13], [14], [15] in which the cut-constructible part of the amplitude is evaluated through 
a particular complex-valued parametrization of the loop momentum. Thanks to the func- 
tional dependence on the complex parameter involved, one is able to separate integral 
contributions to a given cut in combination with OPP subtraction. The rational term is 
evaluated via a loop-level on-shell recursion relation. 

An implementation of the D s -dimensional unitarity cut algorithm [8] (which we will 
call the EGKM algorithm in what follows) in C++ is presented in this paper, very similar 
to that of [16] (written in FORTRAN 95). The goal here is to build up a tool as generic 
as possible that can deal with processes of many external particles in an arbitrary theory 
with fermions, gauge bosons and scalars. When combined with an automated treatment 
of the real radiation (see [17], [18]) this completes the task of evaluating NLO corrections 
to tree-level matrix elements and thus can be used to upgrade to NLO existing matrix 
element generators. The present paper announces the first step in this process, a program 
that calculates numerically loop amplitudes with gluons. 

2. The EGKM approach to loop amplitudes 

The implementation follows the algorithm of [8] which is presented below in short. The 
virtual contribution to QCD amplitudes involving one loop diagrams can be decomposed 
into a sum of terms in each of which color information appears as a factor multiplying a 
gauge invariant quantity called color-ordered amplitude. In the case of N-gluon amplitudes 
the color-ordered amplitudes are amplitudes with a fixed ordering of external legs including 
either a gluon or a fermion loop. In this paper only color-ordered amplitudes with gluon 
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loops are considered. We then have 



[JV/2]+l _ 

A NJu ii= Yl Gr N . c {a)A N , c (2.1) 

C=l <t£Sn /SN:c 

where 

Gr N . c {a) = Tr{t a °m . . . iMc-i) )Tr{t a "^ . . . t a < N ) ) (2.2) 
G riV; i = N c Tr(t a -W . . . £<Miv) ) (2.3) 

and a(i) is some permutation of {1 . . . N}, Sn-,c is the subset of Sn that leaves the ordering 
indicated by the double trace invariant. Given that all A^- c for c > 1 can be written as 
linear combinations of -Ajv,i> we will focus on the calculation of these in what follows. 

A generic color-ordered one loop amplitude with N gluons is just a sum of the relevant 
feynman diagrams. Writing all Feynman diagram contributions under one denominator we 
have 

A *( Pl , P2 ,..., PN ) = j[dl] ^ (2.4) 

where 

di = l 2 di>i = (I + P1+P2 + ■■■+Pi-i) 2 (2.5) 

In the above expression D s denotes the dimensionality of the internal, unobservable, par- 
ticles (gluons in our case). It depends on the regularization scheme and is set to 4 in the 
Four Dimensional Helicity scheme (FDH) [19] or to 4 — 2e in the 't Hooft - Veltman scheme 
(HV) [20]. It is in general different than the dimensionality of the loop momentum (and 
that of the integral), which we denote by D. 

Since the external particles are kept strictly four dimensional, the dependence of the 
numerator on D s is linear for one loop, color-ordered, pure gluonic amplitudes, as can be 
seen by simple numerator algebra: D s enters only as g^ v 9^ v = D s . 

Seen as a function of D H , 



A D S =a° + A 1 -D s (2.6) 

One can, therefore, for a given phase space point and for given external helicities, 
numerically evaluate A Db for two values of D s , determine the D s independent A , A 1 and 
thus get the full A Ds for arbitrary D s . Thereafter, setting D s to 4 recovers the FDH scheme 
and to 4 — 2e the 't Hooft scheme. During all this the dimensionality of the loop momentum 
is kept to arbitrary D with the constraint D < D s . Only after the reduction is performed 
one can set D — > 4 — 2e and evaluate the master integrals as a series in e. 

The problem of numerically evaluating full one loop amplitudes (including the rational 
part) reduces then to the problem of evaluating A Da for two values of D s > D. In the case 
only gluons are present, one can set D s = 5, 6. In the following the calculation of A Ds for 
fixed D s is described. 

It is well known that purely four-dimensional amplitudes can be reconstructed as linear 
combinations of scalar (i.e. with unit numerator) boxes, triangles, bubbles and tadpoles 
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(higher scalar integrals are always reducible to the former [21]): 



Q=\ii,i2,i3Mj T ={11,12,13} 

+ E »»/w5s:+ E (") 

B={ii,i a } 17 11 * 2 S={u} 17 11 

with the summation extending to all possible combinations of propagators (with i,- or- 
dered thanks to the ordering of external legs). Here we use Q, T, B, S as collective indices 
representing ii,i 2 , ■ ■ ■ to keep the formulas readable. 

Equation 2.7 is the result of a reduction process which can be seen as bringing the 
amplitude in the form 

A D.=A( m . _ /*L771 <M0 , /"rjn 



di 1 di 2 di 3 di 4 ... J di 1 di 2 di 3 

^—H,12,13M T=H,12,13 

+ E./ M |f + E>^ 

B=ii,i 2 S=«i 

and then performing the loop integral. 

Note here that the functions (1q(1) ,ct(1) ,&b(0 ,«s(0 depend on I in a particular way. 
For example cIq(1) consists of a piece independent of / and a piece that is proportional to 
the four-dimensional subspace orthogonal to the one spanned by the three momenta that 
enter the propagators. The former is equal to dQ (the coefficient of the scalar integral 
in eq. 2.7) and the latter vanishes upon integration. Terms involving the components of 
I in the subspace spanned by the momenta entering the propagators can be combined in 
propagators and hence appear as constant terms in or in one of ct(0' s - 

When one extends to arbitrary D s > 4 one can also have pentagons 



a^( P1 , P2 , ..., PN ) = Y,[ [di] d . ;fl d . + £/ \diy 



d Q (l) 



H2 ■ j di 1 di2di^di^ 

+E / + E j^W + E / (") 

The equivalent equation for the integrands of the two sides would be 



A°°( P1 , P2 , ..., PN ;l) = J2 , /f] , + e 

+E^ + Eff + E^ (2.10) 
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pen-tuple cuts 

At this point one can multiply both sides of the above equation with five dj's. Then one sets 
the loop momentum such that these five propagators vanish (i.e. multiplies by a suitable 
^-function) and integrates over the loop momentum. 

The left hand side becomes then a product of tree-level amplitudes since the <5-function 
sets the corresponding internal legs on-shell. 

x 5h8e® = / [dl}A D °(pi, P 2, ...,Pn;1) d a d p d^d 5 d e 5 D {1 - I) 

= M^ MX2 M^ MX3 M^ MXi Mf; MX5 Mg^ M (2.11) 

Ai,...,A 5 

with I being the particular D-vector that makes the propagators vanish, and 

M^ MX3 = ^ee(l Xl + P^-upI^pIIW- ■ ■ ^-W-ii-Pk^-i)^)) (2.12) 
where we have explicitly denoted the helicities of external particles, and 

Pk,m = Pk + Pk+l + • • • + Prn (2.13) 

The sums in eq.2.11 extend over all polarization states of the internal, loop propagators 
that have now been put on-shell and have become external legs. Note that the momenta 
I are complex (which is why the three-particle tree-level amplitudes do not vanish). It is 
evident that this step can only be performed if the internal legs have integer dimensionality, 
so D s should be integer and larger than 4. In the case of fermions running in the loop one 
needs also D s to be even. 

Each of the contributing tree-level amplitudes can be evaluated numerically, thus yield- 
ing a value for the left hand side (eq.2.11). 

All terms in the right hand side vanish except the one of the relevant pentagon: 



e E(l) , , , ., , -L> 



di 1 di 2 di 3 di i di 5 



dadpd^dsde 5 (I - I) = e a /3-ySe(l) ( 2 - 14 ) 



so 



^ 7& (l) = e aW£ (f) (2.15) 

If one knows the functional form of e a/ 3 7 5 e (Z) and one has the freedom of finding many 
different Z's that satisfy the unitarity constraints, one can solve for the coefficients of e^Se 
and hence recover it fully. 

Let's call Qi,i = 1...4 the sum of the external momenta between the i'th and the 
i + l'th cut. The five cut propagators can be written as 

d a = l 2 dp = {l+Qif d y = (I+Q1+Q2) 2 d 5 = (I+Q1+Q2+Q3) 2 d e = (I+Q1+Q2+Q3+Q4) 

(2.16) 

Then the solution of the unitarity constraints di = can be found as follows. Find four 
Vi, % = 1..4 such that Qi ■ Vj = dij. This can be done in a straightforward way with the Van 
Neerven-Vermaseren algorithm (see [21], [7]). Then 



= V» + V-^X ( 2 - 17 ) 
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with 

V 1 = Xi vf (2.18) 

and the unit vector in the fifth dimension (i.e. 715 = [0, 0, 0, 0, This automatically 
satisfies d a = as long as n{ = 1. The determined from the other unitarity 

constaints: <%, 7)< 5, e = 0. 

In the pentagon case e(/) = eo is independent of / (see the discussion at [8]), so one I 
value fully determines e Q/ g 7 5 e . 

One should then repeat the process for all pentagon cuts, recovering all eo coefficients. 
At that point, the pentagon contribution to A Ds can be evaluated by multiplying the 
pentagon coefficient with the corresponding scalar pentagon integral. One can do better 
than that, though, since the pentagon is easily reduced in boxes as explained in [8]. Hence 
the pentagon coefficients will be added to the coefficients of the five corresponding boxes. 

Quadruple cuts 

Next, one proceeds to the four-cuts, where a further complication is met: upon taking 
the residue in both sides of eq.2.10 one has, on the right hand side, contributions from 
pentagons that share the same four propagators with the current cut. These contributions 
have to be subtracted from the left hand side to give the particular (1^(1). 

W)=*3k(i)- E f£ fff 1 (2 - 19) 

Let's call Qi,i = 1 ... 3 the sum of the external momenta between the z'th and the 
i + l'th cut. The four cut propagators can be written as 

d a = l 2 dp = {l + Q1) 2 d 1 = {l + Q l + Q 2 ) 2 d s = (l + Q 1 + Q 2 + Q3) 2 (2.20) 

Then the solution of the unitarity constraints d; L = can be found as follows. Find three 
Vi, i = 1..3 such that Qi • Vj = dij and such that n\ ■ Vi = n\ ■ Qi = 0. This vector n\ 
lives in the 4D subspace that is orthogonal to Qi 5 2,3- This can be done in a straightforward 
way with the Van Neerven-Vermaseren algorithm (see [21], [7]). Then 

/> = V^ + ai< + a 5 ng (2.21) 

with 

V» = Xi v? (2.22) 

and the unit vector in the fifth dimension (i.e. 725 = [0, 0, 0, 0, i]). The constraint d a = 
translates into 

a\ + a\ + V 2 = (2.23) 

The other three constraints are used to determine the numerical values of x,- L . 
The functional dependence of the box coefficient on I is, 

d = d + disi + d 2 s 2 e + d 3 s lS 2 e + d^sj (2.24) 
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with 

4 = (*- n *) 2 s i = - n i ( 2 - 25 ) 

i=4..D 3 

being the square of the projection of the loop momentum in the extra dimensions. 

To determine the diS we start with four dimensional loop momenta embedded in a D s 
dimensional space. Then to find do and d\ we need two (four-dimensional) values for I. 

Setting 05 to zero we can get two different Vs: 

Z> 2 = p = V^± \f-\nn p y (2.26) 

from which g?o,i are determined. 
Next, we set 

%\ ai = J-V 1 = a 5 (2.27) 

1%: ai = -^f-V 1 = -a 5 (2.28) 

1% : ai = V-3^ 2 /4 a 5 = ^-V 2 /4 (2.29) 

for the three five-dimensional £3,4,5 that will determine (^2,3,4- Obviously any combination 
of 01,05 with eq.2.23 satisfied would be equally appropriate. 

Triple cuts 

Similarly, for the triple cuts we have 

Calling Qi,i = 1 ... 2 the sum of the external momenta between the i'th and the i + l'th 
cut. The three cut propagators can be written as 

d a = l 2 dp = (l + Qif d 7 = (I + Qi + Q 2 ) 2 (2.31) 

Then the solution of the unitarity constraints di = can be found as follows. Find two 
Vi, i = 1..2 such that Qi ■ Vj = dy and n^, such that • = ■ Qi = 0. Moreover we 
now demand ni • 722 = 0. These vectors ni, ni span the subspace in 4-d that is orthogonal 
to Qi j2 . Then 

/> = + ai< + a2n£ + a 5 n£ (2.32) 

with 

V = XiV? (2.33) 

and rig the unit vector in the fifth dimension (i.e. 715 = [0, 0, 0, 0, i]). The constraint d a = 
translates into 

a? + a\ + + V 2 = (2.34) 
The functional form for c(l) can be chosen to be 

c(l) = C + C1S1 + C 2 S 2 + C 3 (sf - s|) + S1«2(C4 + C5S1 + C 6 S 2 ) + C7«lSe + C 8 S 2 Se + Cgsj (2.35) 
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where 



Si = I • rii (2.36) 

with i = 1, 2 and n\ being the unit vector in the tangent space determined by the momenta 
that enter the triangle vertices. Similarly 

sl= £ {l-nif (2.37) 

i=4..D s 

is the square of the projection of the loop momentum in the extra dimensions. 

To evaluate the first seven coefficients we restrict I to four dimensions. We need seven 
different Z's conveniently chosen, and they are of the form 1 



with 



% = V" + ai, m < + a 2 , m n% 



2ixm 2irm 
a\ = Rcos{ — — ) S2 = Rsin{ — —J m = —3. .3 



(2.38) 
(2.39) 



R = V-V 2 (2.40) 

This setup leads to 

3 

c(U = X ke i2wm/7 m = -3..3 (2.41) 

fe=-3 

where are simple linear combinations of q. The system can be readily inverted to 
yield 

1 3 

A fc = -£e- 2 ™ fc / 7 c(U (2.42) 

-3 

from which Co,. ..,6 are easily obtained. 
Next we set 

al = a 2 = y / -V 2 /2 a5 = y / -V 2 /2 
al = a 2 = -^J-V 2 /2 a5 = y / -V 2 /2 
al = ^-V 2 /2 a 2 = a5 = ^-V 2 /2 (2.43) 

Similarly to the quadruple cut case, any combination that satisfies eq.2.34 is equally apro- 
priate. 

In case \R\ is approaching zero the above definitions become numerically problematic 
and one has to resort to some alternative setup. 

lr The following construction employs a discrete Fourier transform and was employed in [22], [11] as an 
optimized way of constructing and solving the resulting system of equations. 
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Double cuts 

For double cuts we have 

I { i )=X D (i) _ y Weffl _ y (2.44) 

We now have one external momentum Q\ and therefore one v% and three vectors in 
the orthogonal space, n^n^n^. Hence, 

in = yn + ain M + a2n M + a3n M + af . n M (2.45) 

with 

a\ + as + a\ + al = (2.46) 
The functional form of b(l) is (the numerator can contain up to two powers of Z M ) 

b(l) = bo+b 1 s 1 + b 2 S2 + b 3 ss + b 4 (sl-sl) + b 5 (sl-sl) + bQS 1 S2+b 7 s 1 s 3 + bsS2Ss + b 9 s 2 e (2.47) 

with 

Si = l-rii (2.48) 

Restricting ourselves to four-dimensional loop momenta we need nine independent l m 's. It 
is, however, convenient to enlarge the number of equations to ten, so that we can use the 
inverse Fourier transform trick. 

l& = + a 1>m < + a 2im n£ + a 3>ro n£ (2.49) 

with 

„ .2irm. „ ,27rm. 

oi = a 2 = Rcos{— — ) a 3 = Rsin(——) m = -2. .2 (2.50) 
o o 

and 

,27rm. . ,27rm. . . 

ai = Rcos{ — - — J a 2 = a 3 = Rsm( — - — J m = —2. .2 (2-51) 
o o 

with 

R 2 = -V 2 (2.52) 
The last coefficient is evaluated by defining a five dimensional loop momentum as: 

ai = \J-V 2 I2 a 2 = a 3 = a 5 = y / -V 2 /2 (2.53) 

Similarly to the quadruple and triple cut cases, any combination that satisfies eq.2.46 is 
equally apropriate. 

Note that double cuts that isolate a single, external, on-shell line would provide the 
coefficients of massless bubbles with massless, on-shell incoming momentum that vanish in 
dimensional regularization. Hence such cuts need not be considered. 

Single cuts 

Single cuts are not needed in the pure gluonic case (they would multiply gluon tadpole 
scalar integrals which are zero in dimensional regularization). 
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2.1 Master integrals 

Having evaluated all the coefficients of eq. 2.10 for the given phase-space point, we can now 
proceed to perform the loop integration. The constant parts of all coefficients eE,o,dQ t o,CT,OibB,o 
will multiply the corresponding master integrals. The latter are evaluated numerically. The 
sum of all these terms is equal to the cut-constructible part of the loop amplitude. 

a cc = d Q , i Q + ct $ It + bB '° lB ( 2 - 54 ) 

Q T B 

where cIq^ includes <1q : q as well as contributions from reduced pentagons. Any coefficient 
that depends on m-l will vanish upon integration over [dl\. There are, however, coefficients 
that depend only on s^, namely A,cg and 69 of every cut. Those don't vanish upon in- 
tegration. They contribute to the rational part of the loop amplitude. These contributions 
can be reduced to (see section IV. D of [7]) to 

^—E^-E^-E^ <*■«> 

Q T B 

Note that cIq^ doesn't contribute to order in e. 

The scalar integrals needed are all known analytically and have been implemented in 
libraries, first by G. J. van Oldenborgh in FF written in FORTRAN 77 and then by Ellis and 
Zanderighi in QCDLoop written in FORTRAN 90 [23]. For massless theories, there is also 
the implementation [24]. We currently use the QCDLoop library with a wrap code for C++ 
but we plan to interface with the latter as well in the near future. 

3. Summary of the algorithm 

In this section we summarize the reduction algorithm for calculating numerically the virtual 
amplitude for a particular color-ordered set of gluons with fixed helicities in an integer D s 
number of dimensions. One should perform the reduction twice, for two values of D s . In 
the case of gluons D s = 5, 6, but see section 3.2 on a note for D s = 6. 

• Find all the possible pent-uple cuts for the given color-ordered amplitude. 

• for each pent-uple cut 

- evaluate the loop momentum that solves the unitarity constraints (eq.2.17). 

- evaluate the amplitude residue as a product of tree level on-shell amplitudes 
with complex momenta (eq.2.15 - involves assigning polarization vectors to the 
cut legs that carry complex momentum). 

- store the coefficient of the cut. 

- reduce the pentagon to boxes and assign the coefficients accordingly. 

• find all possible quadruple cuts. 

• for each quadruple cut 
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— find the set of loop momenta that solve the unitarity constraints and are conve- 
nient for solving the resulting system (eq.2.26). 

— evaluate the amplitude residue as a product of four tree level amplitudes. 

— subtract the pentagon contributions evaluated at the particular I (eq.2.19). 

— solve the system of equations for do,...,4- 

— store the value of the cut. 

• find all possible triple cuts. 

• for each triple cut 

— find the set of loop momenta that solve the unitarity constraints and are conve- 
nient for solving the resulting system (eq.2.38). 

— evaluate the amplitude residue as a product of three tree level amplitudes. 

— subtract the pentagon and box contributions evaluated at the particular I (eq.2.30). 

— solve the system of equations for co ; ...,9. 

— store the value of the cut. 

• find all possible double cuts. 

• for each double cut 

— find the set of loop momenta that solve the unitarity constraints and are conve- 
nient for solving the resulting system (eq.2.49). 

— evaluate the amplitude residue as a product of two tree level amplitudes. 

— subtract the pentagon, box and triangle contributions evaluated at the particular 
1 (eq.2.44). 

— solve the system of equations for &o,...,9- 

— store the value of the cut. 

• multiply the relevant coefficients with the master integrals to get the final result 
3.1 Polarization vectors 

From the discussion in the previous section it is evident that one needs polarization vectors 
for D s = 5, 6 four- and five-dimensional loop momenta. There are D s — 2 polarization 
vectors in each case. 

When the loop momentum is itself four-dimensional one can use the traditional two 
polarization vectors perpendicular to it and D s — 4 polarization vectors that are unit 
vectors in the extra-dimensional subspace. Care has to be taken, however, since the loop 
momentum is also complex. The only requirement for choosing polarization vectors for 
the internal loop momentum is that they satisfy the proper polarization sums. We employ 
polarization vectors defined as spinor products which have polarization sums corresponding 
to the axial gauge, in that case. 
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When the loop momentum is five-dimensional there are three polarization states. The 
corresponding vectors can be defined in a way analogous to the polarization vectors of a 
massive gauge boson in four dimensions, since then 

l» = l1 D + p^P = il D -l> = Q (3.1) 

The polarization vectors corresponding to 1% D with mass I 2 are the same as those corre- 
sponding to in 5 dimensions. 

3.2 Tree-level amplitudes 

Since one can avoid entirely a fully six-dimensional loop momentum (see the discussion 
on the necessary loop momenta of section 2), it is easy to see that the six-dimensional 
tree-level amplitudes needed for the construction of the cut residues always decouple into 
five-dimensional tree- level amplitudes (when the polarization of the loop legs are ±,5) or 
amplitudes with scalars in the place of the loop legs (when the loop legs are both polarized 
in the 6-th dimension). The tree- level amplitude can always be written as 

A 6 ( Pl ,e 1 ;...; Pn ,e n ) = 4M%e v n (3.2) 

Since M^ v is purely (at most) five-dimensional 2 , if e\ and e n are polarized along the 6th 
dimension (i.e. they are equal to [0, 0, 0, 0, 0, 1]) the only non-zero term will occur for 

K> "> A6 sc9»» (3-3) 

where A Q SC is the part of the subamplitude proportional to g^. Therefore one really needs 
tree-level amplitudes with gluons in D s = 5 and tree-level amplitudes with gluons and 
scalars in D s = 5. 

Tree-level color-ordered amplitudes with gluons can be calculated recursively using 
the Berends-Giele recursion relation [25]. One needs look-up tables to avoid calculating 
the same subamplitude more than once. Alternatively, one can use a bottom-up Berends- 
Giele approach, as employed in [26] in which it is by construction guaranteed than no 
subamplitude is calculated more than once. We have experimented with the bottom-up 
approach with explicit code for each number of external legs generated and compiled once. 
It is slightly faster for N < 12 but becomes slower for N > 12 presumably due to code 
bloating. Here we use the top-down apporach. The gluon-scalar amplitudes are similarly 
generated. 

Since the five-dimensional gluon and scalar color-ordered amplitudes are the building 
blocks of this algorithm, efficient evaluation of them is of paramount importance. In table 1 
we give the CPU time needed to evaluate a given tree-level color-ordered amplitude in the 
gluon and scalar case for a number of external legs varying from four to twentytwo. Typical 
cpu times needed for the calculation are estimated by averaging over the time needed for 
the evaluation of 10 5 phase space points generated by RAMB0 [27] with a given set of cuts 
(\m\ < 3,p T ,i > QMy/s,Rij = J(j>lj + rifj > 0.4) on an Intel Xeon X5450 @3.00GHz. 



2 It is five-dimensional when the loop momentum is five dimensional, and four-dimensional otherwise. 
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Time for a tree amplitude 


N 


t(fis) 


N 




N 


t(/xs) 


4 


5 


11 


178 


18 


1326 


5 


9 


12 


252 


19 


1649 


6 


17 


13 


351 


20 


2032 


7 


29 


14 


475 


21 


2482 


8 


50 


15 


631 


22 


3004 


9 


79 


16 


818 






10 


121 


17 


1048 







Table 1: CPU time averages for tree-level amplitudes over 10 5 phase space points generated by 
RAMB0 [27] with a given set of cuts (|nj| < 3,p T ,j > 0.01y/s,Rij = \J~$\j + «|j > 0.4) on an Intel 
Xeon X5450 03.00GHz. 



4. Numerical issues 

Issues of numerical stability have to be addressed in any algorithm performing numerical 
integration of one-loop amplitudes, due to the possible appearance of vanishing (or almost 
vanishing) Gram determinants in one or the other part of the calculation. In this particular 
algorithm the problem can appear in the evaluation of the reduction coefficients, and, 
in particular, when calculating the solution to the unitarity constraints for a given cut. 
Moreover, the subtraction of the contribution of higher cuts in quadruple or lower cuts is 
another potential source of loss of accuracy. 

The way to deal with this problem employed here is to check for any given point, 
the coefficient of the poles in the e-expansion, known a priori in an analytic form thanks 
to [28], [29], [30]. If the (normalized) difference between the calculated pole coefficient and 
the analytically evaluated one is more than a predetermined value (say 10 ) the phase- 
space point is considered problematic. Thereafter the whole algorithm is performed again 
using quadruple precision complex numbers (implemented in the QD library [31]). The 
quadruple version of the code is of course much slower than the native double precision 
version, so one hopes that the fraction of times it needs to be employed is relatively small. 
In the case of six gluons approximately 5% of the points need quadruple precision if the 
switch accuracy is set to 10 -4 . 

In figures 1,2 we show the log of the relative precision, 

dE 1 2 = — - Eanalytic ^ 2 (4.1) 

Eanalytic,l,2 

for the double and single poles for the double precision version and for the version with the 
quadruple precision switched on. A similar distribution for the relative error is observed 
when the number of gluons increases, albeit the peak of the distribution shifts towards 
larger errors. This is expected to pose serious problems when N > 12, but given the 
fact that one can even resort to arbitrary precision arithmetics (APREC [32]), cases of 
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phenomenological importance should always be manageable with a relatively small penalty 
in time. 

Note also that, despite the fact that relative errors for the poles are correlated with 
each other and with the relative error for the finite part, the latter is not guaranteed to 
behave well when the former do, due to potential instabilities in coefficients of finite scalar 
integrals. There is a variety of checks one can perform to try to detect such cases so that 
quadruple precision arithmetics is switched on. 



5. Performance 




1 1 1 1 1 1 1 1 1 1 1 

4 6 8 10 12 14 16 18 20 22 

Number of gluons 

Figure 3: CPU time for the color-ordered, helicity- fixed, N gluon amplitude in double precision 
on an Intel Xeon X5450 @3 . 00GHz 
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5 


12 


10 


1193 


15 


32465 


20 


564330 


6 


34 


11 


2376 


16 


59145 


21 


1011710 


7 


91 


12 


4809 


17 


107145 


22 


1848280 


8 


222 


13 


9380 


18 


190635 







Table 2: Typical cpu times needed for the calculation arc estimated by averaging over the time 
needed for the evaluation of 10 5 phase space points generated by RAMB0 [27] with a given set of cuts 
(\rii\ < 3,p T:J > 0.01y/s,Ri tj = yjtplj + nfj > 0.4) on an Intel Xcon X5450 03.00GHz. 

The code is presently able to calculate the full virtual part of color-ordered amplitudes 
with gluons, up to an arbitrary number of external legs and for arbitrary helicities, including 
the 1/e 2 and 1/e coefficients, as well as the rational part of the finite term in the e- 
expansion. 

Checks with published results for specific phase-space points at [16] and [8] have been 
performed successfully. 

Typical cpu times needed for the calculation are estimated by averaging over the time 
needed for the evaluation of 10 5 phase space points generated by RAMBO [27] with a given 
set of cuts (|n;| < 3,pr,i > 0.01y/s,Rij = \J 4>\j + n \j > 0.4) on an Intel Xeon X5450 
@3.00GHz. Table 2 shows those averages for N = 4. .22. In figure 3 these averages are 
plotted together with the equivalent tree- level ones multiplied by 10 3 . 

6. Outlook 

We have presented here the first step in an effort to produce a generic implementation 
of the EGKM algorithm [8] for calculating numerically one-loop amplitudes. The current 
version involves gluons as external and as virtual particles. The necessary building blocks 
are color-ordered multi-gluon tree-level amplitudes and tree-level amplitudes with gluons 
and scalars in five dimensions. They are all evaluated via a Berends-Giele [25] recursion 
relation. 

The performance of the code is quite satisfactory and agreement has been reached with 
previously published results. Numerical stability issues are all addressed by resorting to a 
quadruple precision version of the code. 

There are ample phenomenological applications for a generic code that includes fermions 
and electroweak gauge bosons, all within the immediate reach of the D s -dimensional uni- 
tarity cut approach to OPP reduction, as demonstrated explicitly in [10]. Extending the 
implementation to (massless or massive) fermions and gauge bosons is relatively straight- 
forward. To this generalization we turn in the future. 
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